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Resonant alignment of microswimmer trajectories in oscillatory shear flows 
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Oscillatory flows are commonly experienced by swimming microorganisms in the environment, 
industrial applications and rheological investigations. We experimentally characterise the response 
of the alga Dunaliella salina to oscillatory shear flows, and report the surprising discovery that algal 
swimming trajectories orient perpendicular to the flow-shear plane. The ordering has the charac¬ 
teristics of a resonance in the driving parameter space. The behaviour is qualitatively reproduced 
by a simple model and simulations accounting for helical swimming, providing the mechanism for 
ordering and criteria for the resonant amplitude and frequency. The implications of this work for 
active oscillatory rheology and industrial algal processing are discussed. 

PACS numbers: Valid PACS appear here 


Many swimming microorganisms experience shear flow 
in natural and industrial processes. Swimming is strongly 
biased by environmental cues and fluid shear [2, 2|, with 
significant implications for ecology [I] and industrial ex¬ 
ploitation [ 4 ]. Classic examples include directed swim¬ 
ming relative to light (phototaxis) and hydrodynamic fo¬ 
cusing in down-welling flow due to viscous and gravita¬ 
tional torques (gyrotaxis) HI- 

There is great potential to use individual and collective 
microswimmer behaviour to improve microbial biotech¬ 
nology, such as algal photobioreactor design [dj. For 
example, gyrotactic microorganisms in laminar channel 
flow tend to focus and so drift faster and diffuse less 
than non-swimming cells or nutrients Hi while in tur¬ 
bulent flows cells accumulate in transiently downwclling 
Si I or strongly accelerated 0 regions. Horizontal 
shear flows can trap gyrotactic swimmers (a mechanism 
for oceanic thin layers) Q and modify hydrodynamic in- 
stabilites and patterns (bioconvection) [11, |l2|. Photo¬ 
taxis and shear flow can combine to drive cell focusing 
[bl and pattern formation 0- Complex transport dy¬ 
namics can even result from relatively simple shear flow 
00. The rheology of active media is also of recent 
interest: suspensions of swimming bacteria behave less 
viscously jl6l j and algae more viscously 0 than dead 
cells. 

Here, we investigate the interaction of the swim¬ 
ming alga Dunaliella salina with oscillatory shear flows. 
Surprisingly, in experiments swimming trajectories are 
strongly ordered by the flow for particular driving param¬ 
eter values. The ordering mechanism is distinct from that 
observed recently with Dunaliella primolecta, with con¬ 
stant, strong shear flows |lSj. We explore the observed 
resonant ordering employing simple but predictive mod¬ 
els, and discuss implications for active oscillatory rheol¬ 



FIG. 1. (a) Schematic of apparatus to image suspensions 

of algae in oscillatory shear flows, (b) Flow oscillation axis 
e y (not always aligned with the image x axis, but consistent 
across repeat experiments) and the angle 9 of the swimming 
trajectory displacement Ar to the x-axis. 


ogy and industrial processing of swimming algae. 

D. salina (CCAP 19/18) cells were grown on modi¬ 
fied Pick medium 0,j2(j] under 12:12 light/dark cycle at 
21°C. All experiments were carried out at this temper¬ 
ature. Cells were concentrated by upswimming (gravi- 
taxis) using cotton wool rafts 0 . Dilute (10 6 cells/ml) 
suspensions were subjected to oscillatory shear on the 
stage of an optical microscope (Olympus BX51). The 
suspension was placed between two transparent parallel 
plates 400 /im apart, the top plate connected to an elec¬ 
tromechanical drive that sinusoidally sheared the suspen¬ 
sion (Figure [1]). Plate parallelism was ensured by zeroing 
sample capillary flow. Video sequences of sheared algae 
were acquired using a Mikrotron MC1310 at 10 x (NA 
0.25) using red-filtered bright-held illumination to min¬ 
imise phototaxis [2lJ . Sequences were captured in a plane 
equidistant from top and bottom plates at depth 200 /im. 
Algae were tracked using MATLAB versions of estab¬ 
lished algorithms [22|, [23J. The direction of the imposed 
oscillating flow was inferred from short-time cell trajec¬ 
tories and confirmed by tracking PEGylated polystyrene 
colloids (Supplementary Materials, FigSl). Cell observa- 
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tion before and after measurements found that the appa¬ 
ratus did not damage the cells (e.g. deflagellation). 

With no flow, D. salina swimming trajectories were 
distributed isotropically in the horizontal plane (Figure 


HO- As in [24| , gravitactic bias was not evident for tracks 


in this plane on experimental timescales. On application 
of oscillatory shear with amplitude A and frequency / 
trajectories might be expected simply to reflect superpo¬ 
sition of isotropic swimming and oscillatory advection. 
This was indeed the case for some driving parameters, 
such as A — 120 (im and / = 6 Hz (FigOi)- However, 
for A = 224 /im and f = 2 Hz swimming trajectories un¬ 
expectedly aligned perpendicular to the flow-shear plane 
(Fig[2t)- Alignment can be quantified by the start-to- 
end displacement vector, Ar(r), of the trajectory of each 
swimmer (see[l]o). Ar(r) evaluated at short times within 
oscillation cycles provides the oscillating flow direction 
e y ~ Ar(T —>• 0). However, evaluating Ar(r) over the 
largest available integer multiple n of the oscillation pe¬ 
riod, r = n/f, excludes cycle-by-cycle flow oscillations 
and provides the orientation 9 of the swimming trajec¬ 
tory; i.e. it measures net swimming progress. 

Distributions P{&) of orientations are presented beside 
the trajectories in Figure [2] The distribution for / = 2 
Hz, A = 224 (im shows how trajectories orient along 
the line perpendicular to the flow-shear plane, but are 
equally likely in either direction along this line. If the 
flow is halted, the distribution returns to uniform (Figure 
SI, Supplementary Materials). Swimming speed distri¬ 
butions, also shown in Figured appear to accentuate the 
importance of strong swimmers during alignment. How¬ 
ever, for all cases modal speed is «60 pm s _1 , which 
compares well with recent measurements Q. 

To statistically quantify observed alignment as a func¬ 
tion of the flow parameters, we count trajectories with 
displacement Ar oriented perpendicular, N±, and paral¬ 
lel, AT||, to the flow direction and define R = Nj_/N\\. By 
parallel (perpendicular) displacements we mean swim¬ 
ming orientations within ±7r/4 of the flow (vorticity) 
axis. The surface plot in Figure [3] illustrates how the 
alignment has the characteristics of resonance, occupy¬ 
ing a small region of the imposed flow parameter space. 
Fixed amplitude and frequency sections of the ordering 
surface are shown in Figure [Hand discussed below. 

Discussion. The observed trajectories result from the 
combination of shear flow and swimming. Many algae 
swim helically, in part to facilitate phototaxis via a di¬ 
rectional eyespot [26|. Thus we model the algae as he¬ 
lical swimmers in a flow. Following [l|, we assume that 
a cell at position r swimming with speed v in direction 
p has an intrinsic angular velocity u>h n about an axis 
n, where p n = cos/3, for constant angle /3, due to 
an asymmetric (non-planar) flagellar stroke. Hence, if 
the only external torque is due to flow with velocity u 
and vorticity u>, spherical swimmers obey r = u(f) + vp, 
p = (icu(f) + u>h n ) x P an d n = x n. Experimen- 
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FIG. 2. Trajectories, orientation distributions and speed of 
D. salina. (a) No flow: isotropic trajectories, (b) Oscillatory 
shear flow with amplitude A = 120 pm and frequency f — 6 
Hz: tracks show oscillation at short times, but retain isotropy 
of swimming orientation over longer times, (c) A = 224 pm 
and / = 2 Hz: swimming directions align perpendicular (blue) 
to the flow oscillation direction (cyan). Pink/green circles 
denote start/end points of each 5s track. 


tal observations suggest that the fluid velocity is u (t) = 
7 oo-Z’ cos(uidt)e x , with vorticity u>(t) = 7 oo cos (ujdt)e y . 
Here e x and e y are along the positive flow and vorticity 
directions, respectively, 700 = u>dA/H is the maximum 
shear rate and u>d = 27r/ is the angular driving frequency 
(recall A. / and H are driving amplitude, frequency and 
gap width, respectively). Nondimensionalising lengths 
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such that p = n = (sin0sin<I>,cos0,sin0cos<I>), where 
0 increases from the direction of vorticity, e y , along the 
T-axis, and $ is measured from the Z -axis [e.g. (0, <f>)= 
(tt/2, 7 t/2 ) is along the X-axis]. Equations (IA.131 l2|) give 


X = sin $ + F cos {t)Z (3) 

y = »\\ ( 4 ) 

Z = v l cos $ (5) 

0 = 0; 4> = — cos (t) (6) 


where u± = v sin 0o and = v cos 0o are the nondimen- 
sional swimming speed components perpendicular and 
parallel to e y . Integration of (1A. 171) and (IA.19I) yields 
Y(t) = Yq + zzyi, 0(f) = 0 O and 


FIG. 3. Alignment ratio R as a function of driving amplitude 
A and frequency / for oscillatory flow in a shear cell with H = 
400 fi m gap width. A is a measure of horizontal displacement. 


with A and times with 1/uid, the model equations read 
r = T cos (t)Ze x + v p (1) 


P = 


■cos(t)e y +Q i n 


x p; n 


■ cos (t)e y x n (2) 


where T = , y 00 /ojd = A/H is the dimensionless shear 
rate/amplitude, = Ud/ijJh is a frequency ratio and 
v = v/(Au>d) = iV(rn) a nondimensional swimming 
speed, with vq = v/(uihH). The Cartesian representa¬ 
tion of (IA.13P ) , with Z measured from the bottom plate 
and XY the flow-vorticity plane (figure©, was solved 
numerically to simulate a suspension of swimmers (see 
Supplementary Materials). The ratio R was computed 
from many simulations with a uniform distribution of 
initial swimmer orientations. 

Non-swimmers (zz —>• 0) passively follow the imposed 
oscillatory flow. In the absence of flow, the governing 
equations predict helical trajectories: p rotates around 
n (the helix axis) with frequency u>h- With the experi¬ 
mental resonance close to helical swimming frequency, it 
is tempting to think that helical trajectories are respon¬ 
sible for the observed alignment. This is only partially 
true, as we shall see. 

First we ask if the combination of ‘non-helical’ swim¬ 
ming and oscillatory shear alone is sufficient to induce the 
observed behaviour. Figure [4] displays the ordering ra¬ 
tio (experiment a,b; simulation c,d) and simulated swim¬ 
mer trajectories (e) for a model where helical swimming 
is switched off (dashed lines). While R peaks at char¬ 
acteristic values of the driving amplitude T (Fig [4jc), it 
is entirely independent of driving frequency (Fig [4]:1) in 
stark contrast to the experiments. The non-helical limit 
does, however, predict alignment: it is worth considering 
further. In this limit (/? —>- 0), the governing equations 
simplify considerably if we choose Euler angles 0 and $ 


$(f) = $o + ^sin(t). (7) 

(Recall that T = A/H is the non-dimensional shear rate.) 
As the Y component of the trajectory grows linearly in 
time, independent of shear, alignment can only depend 
on the coupled X and Z dynamics. In particular, closed 
orbits in the Xzf-plane are present at resonance, see fig- 
ure©, panel (i). For such orbits, progress only in the In¬ 
direction is possible, leading to alignment. Off-resonance, 
orbits are open and cells can progress in X and Z direc¬ 
tions (figure©, panel (ii)). 

This phenomenology can be understood in terms of 
oscillatory Jeffery dynamics of the swimmer orientation. 
We see from © that oscillatory shear forces swimmer 
orientation in the vertical X,Z-plane to describe circular 
arcs swept sinusoidally in time (contrast this with circu¬ 
lar Jeffery orbits in steady shear flow, uid —> oo), with 
angular amplitude T/2. Folded orbits only arise when 
shear is sufficiently large to rotate swimmer orientation 
by integral multiples of 7r, so it can make no net progress 
during a cycle (see©:). This provides a prediction for the 
resonant ordering amplitude T res ~ 2irn, n £ Z, in good 
agreement with non-helical simulations (Fig ©). The 
latter agree qualitatively with the experimental results 
in Fig ©. Quantitatively, smaller values of T are suffi¬ 
cient to induce ordering in experiment. This could be 
due to second order mechanisms neglected in our simple 
model, as discussed below. 

The model for non-helical swimmers predicts align¬ 
ment as a function of amplitude, but does not reproduce 
the experimentally observed dependence on both driving 
amplitude and frequency (FigHr). D. salina is a helical 
swimmer, rotating at 1.5 — 2 Hz 4|. This second, internal 
frequency provides the possibility of further resonance. 
Indeed, with /3 ^ 0 (recall j3 is the angle between p and 
n) numerical results reveal that ordering is frequency de¬ 
pendent (Figure©!). Position equations are unchanged, 
but depend on more complex orientation dynamics re¬ 
sulting from the coupling of flow-induced and helical ro¬ 
tation. Cell orientation angles 0 P and defining p, 
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evolve according to (see Supplementary Materials) 

0 p = fl _1 sin0 n sin($„ - $ p ), (8) 

= 77 cos (t) (9) 

+ fi -1 [cos0„ — sin0„ cot 0 p sin(4>„ + $ p )], 

whereas angles for n satisfy (IA.19|) . such that 0„ = 0, 
= r cos(t)/2. As the p-dynamics are slaved to the n- 
dynamics, trajectories with helical swimming do retain 
broad features of non-helical orbits in oscillatory shear, 
see Figure®, but they are nevertheless qualitatively per¬ 
turbed (even for infinitesimal j3). Thus only for particu¬ 
lar driving frequencies and amplitudes does helical swim¬ 
ming produce alignment-inducing orbits. The frequency 
condition for ordering can be obtained by considering the 
case of a swimmer with n in the direction of vorticity: 
0„(O) = 0 = $„(0); 0p(O) = /?; and <f> p (0) = 0. Equa¬ 
tion (IA.21I) then integrates to <f> p (f) = ^ sin(t) + a 

helical phase can perturb simple Jeffery rotation by the 
flow. Only when Q ~ 1, i.e. when the driving and he¬ 
lical phase are synchronised, can a resonant value of T 
give alignment, agreeing with both simulations and ex¬ 
periments (Fig®>,d). 

Oscillatory Jeffery orbits and helical swimming provide 
a first order explanation for the dependence of resonance 
in experiments on both frequency and amplitude. The 
model is predictive: swimmers with different helical fre¬ 
quencies should display a different resonance spectrum 
(Fig lUi, inset). This description could be extended to in¬ 
vestigate how additional effects such as taxes, orientation 
noise, inertia and cell shape either dominate, compete or 
act in concert with helical swimming to affect resonance. 
Intriguing possibilities include stochastic resonance due 
to noise in the flow velocity gradients [§8|, the active 
response of flagella to shear [2 0] , and effective shape 
and response to shear due to flagellar beats [3fj. A full 
analytical investigation of the nonlinear dynamics of the 
helical model in oscillatory shear is beyond the scope of 
this paper, but it is clear that much is to be discovered, 
analogous to structures observed for swimming cells in 
Poiseuille flow jl5j ]. 

Conclusions. We have demonstrated the surprising 
response of swimming microalgae to oscillatory shear 
flows, producing an alignment of trajectories with a set 
of resonance peaks in the parameter space of driving fre¬ 
quency and amplitude. A simple model combining shear 
and non-helical swimming predicts resonant alignment 
of the trajectories of swimming cells, but only when he¬ 
lical swimming is included does the experimentally ob¬ 
served frequency dependence of the resonance arise. The 
rich dynamics and the counterintuitive interactions be¬ 
tween swimmers and flow revealed by our experiments 
and modelling have implications for both active suspen¬ 
sion rheology and the design of novel cell processing 
methods. While simplified models of swimmers (rod or 
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FIG. 4. Alignment ratio A as a function of amplitude Y 
and frequency Yl from experiment (a, b) and simulation 
(c, d). Simulation predictions are shown with (/3 = tv/6, 
solid) and without (/3 = 0, dashed) helical swimming (fre¬ 
quency fh = 2Hz), only the latter providing frequency de¬ 
pendence. Trajectories are also shown in (e): at resonance, 
[T,fi] = [7.5,0.75], for XZ plane (i) and 3D views (ii); off- 
resonance [r,fi] = [4,0.75] (iii). Only trajectories closed in 
XZ provide ordering. Inset of (d): simulations of Pfiesteria 
piscida in oscillatory flow. When predating, this dinoflagel- 
late increases its helical frequency from 1 (solid line) to 2 Hz 
(dash dot) [3lj. We predict a measureable shift in the reso¬ 
nance peak. 


spheroidal pushers and pullers with no helical motion) 
appear adequate to explain active rheological phenom¬ 
ena such as shear-thickcning in algal suspensions (T?! [32j , 
the current work suggests such models may fail in ac¬ 
tive oscillatory rheology experiments. Biotechnologically, 
the results hint at methods for im pro vement in efficiency 
of the algal processing pipeline [331 ]. For example, in 
downstream processing of useful microalgae, like the /3- 
carotene producer D. salina , cells commonly experience 
oscillatory and squeezing flows [.'i ll . Resonant alignment 
will provide boundary accumulation over times L/v, 
where L is the size of the shear plate and v is the swim¬ 
ming speed. This may be ‘engineered out’ by tuning pro¬ 
cess parameters from resonance; or it may be fruitfully 
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exploited as a new way to guide and harvest cells. 
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SUPPLEMENTARY MATERIALS 
Shear effects on swimming and flow 

In order to establish that the swimming of Dunaliella 
salina in the absence of flow was not affected by exposure 
to oscillatory shear, the trajectory orientation and swim¬ 
ming statistics of cells were monitored before and after 
every shearing experiment, see Figure 0). Isotropy and 
mean speed return to normal once the oscillatory flow is 
halted. The only effect is a drop in the total number of 
swimmers in the field of view, probably a result of cell 
loss in the field of view generated by the anisotropic mi¬ 
gration caused by the oscillatory shear (as discussed in 
the conclusions of the main text). The same 0 shows 
the simultaneous tracking of swimming algae and PEGy- 
lates polystyrene colloids, verifying they passively follow 
the flow, as expected. 


Simulation model 

As mentioned in the main text, mechanisms for the 
experimental trajectory ordering of Dunaliella salina in 
oscillatory shear flows were explored by numerical simu¬ 
lations of helical swimmers. The equations of motion of 
a swimmer were based on a recently proposed model by 
Bearon [lj. In this model (notation as in the main text), 
cell position r evolves by swimming with mean speed v 
along the direction p. Through an intrinsic torque due 
to nonplanar swimming the cell also generates rotation 
about a direction n. The dynamics of swimmers in a flow 
u is governed by 


r = u + up 

(A.10) 

p = w c xp 

(A.ll) 

n = uj c x n 

(A.12) 


where the angle between p and n is fixed and given by 
p n = cos (3 and uj c is the net angular velocity of a 
cell resulting from the balance of torques on it. Bearon 
considered external torques on the cells due to gravity 
on bottom heavy spherical swimmers in a flow, so that 
u> c = ujgP x e~ T + 1 u;, where u> g is the gravita¬ 
tional reorientation frequency, u>h is the intrinsic angu¬ 
lar speed of helical rotation and ui is the flow vorticity 
(horizontal and vertical linear shear flows were solved 
for in El)- The model neglects known effects, such as 
cell shape and stochastic reorientation. These could eas¬ 
ily included (e.g. as in 1 ), but to capture a minimal 
mechanism for our observations we neglect them here, 
and additionally neglect gravitational torques. These 
torques will only act on cells swimming with a signif¬ 
icant vertical component. For these, gravitational re¬ 
orientation of the cells acts on timescales l/ui g ~ 10 s 
[3], which are large compared to the flow timescales for 


most driving frequencies considered. Thus, flow reorien¬ 
tation is expected to dominate the dynamics, as verified 
from simulations (see Figure 0. We then assume the 
simple oscillatory flow u(f) = ^^Z cos{u)dt)e x (vorticity 
u>(t ) = 7 oo cos (u!dt)e y ), where 700 = aJdA/H is the maxi¬ 
mum shear rate for driving angular frequency ujd = 27T fd 
and amplitude A. Nondimensionalising in terms of these 
time (1 /ujd) and length (A) scales, equations (IA.10HA.f2l) 
reduce to equations ( 1 ) and ( 2 ) of the main text 


r = F cos (t)Ze x + up 


P = 


— cos(f)e y + 1 n 


x P 


n = — cos (t)e v x n 


(A.13) 
(A.14) 

(A.15) 


where, as in the main text, we define the dimensionless 
driving amplitude/shear rate F = A/H and frequency 
U = ujd/tjJh- The parameter v = v/(Au>d ) gauges the rel¬ 
ative magnitude of advection by swimming to advection 
by flow. Re-writing v = r'o/fir, with uq = v/(Hujh) a 
constant, we see that v cannot be varied independently 
of the dimensionless driving parameters. In [lj, where 
gravity can play a significant role, the ODEs for the com¬ 
ponents of (IA.1311A.T5l) were derived in a Cartesian coor¬ 
dinate frame where Euler angles are defined with respect 
to the vertical. We consider instead Cartesian frame with 
Euler angles 0 and $ describing the helical rotation di¬ 
rection defined from the vorticity direction e y along the 
Y-axis: n = (sin 0 n sin $ n , cos 0 n , sin 0„ cos $„), as de¬ 
scribed in the main text. This choice simplifies the angu¬ 
lar equations of motions and makes the dynamics more 
physically transparent.The equations of motion for the 
swimmers then become 


X = ^sin0„ sin'll + Fcos(wdf)Y (A.16) 

Y = VCOS&n (A.17) 

Z = v sin 0„ cos (A.18) 

0„ = 0; <F„ = - cos(f) (A.19) 

0 p = fl -1 sin 0„ sin(d> ri — <F p ) (A.20) 

= F cos(f )/2 


+ 1 [cos 0 „ — sin 0 „ cot 0 P sin($ ra + $p()jk.21) 

In the limit /? —>• 0, helical effects are negligible. In 
this limit, losing subscripts and defining u± = ^sin 0 o, 
v\\ = v cos 0o, equations (1A.161 [/Oil) reduce to equations 
(3 — 6 ) in the main text, where the dynamics for 0, $ and 
Y are trivially solved for analytically. The non-integrable 
equations of system (IA.16I IA.21I) for swimming cells were 
solved numerically with MATLAB (Mathworks, Nat¬ 
ick, MA, USA) using a Runge-Kutta-Fehlberg (RK45) 
method in a 1800 x 1800 x 100/nn periodic box. Realistic 
D. salina mean swimming parameters (along-helix speed 
v = 80/rm/s and intrinsic frequency fd = 2Hz) were 
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(c) 


FIG. 5. Trajectory orientation and speed PDFs before (a) and after (b) oscillatory flow (A = 224 /im, / = 2 Hz as in figure 2c 
of main text). When oscillatory shear ceases, trajectory isotropy is recovered, (c) Swimmer trajectories (colour) and colloids 
(black) for A = 150 /.tm, / = 2 Hz. The latter passively follow the oscillatory flow. 
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FIG. 6. Alignment ratio simulations with and without gravi- 
taxis, which makes little difference to ordering dominated by 
flow and helical contributions. 

used Q, ignoring phenotypic variations across the popu¬ 
lation. The helical angle between the swimming direction 


p and the helical axis of rotation n can be estimated from 
[3 = arccos(u/u p ) = 0.5, where v p = 60/im/s is the mean 
progressive speed along the helical axis [4j. Cells were 
set off at Z = H/ 2, with initial values of p 2D-isotropic 
in the horizontal plane and n to be rotated by (3 out of 
the plane. Realistic driving parameters (fa = 1 — GHz 
and A = 50 — 3500/im), and a fixed gap width 100/im 
were used. The ordering ratios R were obtained from the 
simulation position data, as for experiments. 
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